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Abstract 

Hydrodynamic equations are used to identify the final state readied by a freely evolving granular 
gas above but close to its shear instability. The theory predicts the formation of a two bands shear 
state with a steady density profile. There is a modulation between temperature and density profiles 
as a consequence of the energy balance, the density fluctuations remaining small, without producing 
clustering. Moreover, the time dependence of the velocity field can be scaled out with the squared 
root of the average temperature of the system. The latter follows the Haff law, but with an effective 
cooling rate that is smaller than that of the free homogeneous state. The theoretical predictions 
are compared with numerical results for inelastic hard disks obtained by using the direct Monte 
Carlo simulation method, and a good agreement is obtained for low inelasticity. 

PACS numbers: 45.70.Mg, 45.70.Qj, 47.20.Ky 
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I. INTRODUCTION 



Granular gases are ensembles of macroscopic particles whose dynamics is controlled by 
inelastic binary collisions, separated by ballistic motion The forces are usually short 
range and repulsive and, therefore, in the simplest form granular gases are modeled as a 
collection of smooth inelastic hard spheres or disks. In recent years, the study of granular 
gases has attracted a lot of attention , among other reasons because they serve as a 
sensitive proving ground for kinetic theory and non-equilibrium statistical mechanics. 

A peculiarity of granular gases is their tendency to spontaneously develop pattern forming 
instabilities when evolving freely. This includes the so-called shearing and clustering insta- 
bilities 0, [5I, in which the system develops patterns both in the flow field and the density 
field. Analysis of the hydrodynamic Navier-Stokes equations for granular gases indicates 
the existence of several possible scenarios explaining the development of these instabilities 
and their mutual influence. An illuminating discussion of them is given in ref. ^. For 
instance, the existence of non-stationary channel flows in freely cooling gases leading to an 
attempted finite-time density blowup, which is arrested by heat diffusion, has been shown 
recently 0, ^. 

In this paper, a shear flow state of freely evolving dilute granular gases will be described. 
Its origin is tied in with the shear mode instability and the analysis to be presented here is 
restricted to systems whose size is slightly larger than the critical size for this instability. The 
latter happens to be smaller than the critical size for the instability of one of the longitudinal 
modes associated with the density [9|, HOj . Consequently, the spatial inhomogeneities of the 
hydrodynamic fields are expected to be small. Although the flow is non-stationary, all 
the time dependence can be expressed through the spatially averaged temperature and, 
therefore, it can be described in terms of dimensionless time-independent quantities. More 
precisely, it is an inhomogeneous two band shear state. It must be stressed that this state 
is exhibited by a system with, for instance, periodic or elastic boundary conditions, without 
bein g d rive n by the boundaries as it happens with the so-called simple or uniform shear 
flow jlll. I12I. Il3[|. Non-stationary states in which all the time dependence occurs through the 
temperature, like the homogenous cooling state IJ], as well as non-uniform states where all 
the spatial dependence also occurs through the temperature are peculiar of granular gases 

Q. 
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The existence of this free shear state has been suggested previously in the framework of a 
time-dependent Ginzburg-Landau model for granular gases [l^, and also as the result of a 



nonlinear analysis of the shearing instability 



16| . The theory presented in this paper differs 



from those studies both in the method and the results, where rather strong discrepancies 
occur as it will be discussed along the paper. The interest here is focussed on the complete 
identification of the hydrodynamic fields characterizing the free shear state and the compar- 
ison of the theoretical predictions with the results of numerical simulations of the particles 
composing the system. Of course, this provides a very demanding test of the validity of the 
macroscopic description provided by the hydrodynamic Navier-Stokes equations for granular 
gases. 

The remainder of the paper is organized as follows. In Sec. [ITl the nonlinear Navier-Stokes 
hydrodynamic equations for a dilute granular gas are shortly reviewed, and particularized 
for the free shear state. This is macroscopically characterized by presenting gradients only in 
one direction and by a velocity field perpendicular to the gradients. Moreover, it is assumed 
that the density inhomogeneities are small. By introducing appropriate dimensionless length 
and time scales, a closed equation for the velocity flow is derived in Sec. IIIII The solution of 
this equation shows that the amplitude of the velocity field decays monotonically in time. 
Moreover, by analyzing the temperature equation, an explicit expression for the amplitude 
of the steady density fluctuation is obtained. Also, it is seen that the ratio between the 
average temperature of the system and the square of the amplitude of the velocity field is 
time independent. 

In Sec. IIVI the simulation method used is indicated. It is restricted to low density 
gases, whose one-particle distribution function obeys the Boltzamnn equation, but it must 
be stressed that it is an A^-particle simulation method, therefore providing information on 
all the space and time correlations. It is observed that the system actually tends to a 
steady situation having the assumed properties of the free shear state. Then, the measured 
hydrodynamic quantities are compared with the analytical theoretical expressions derived 
from the Navier-Stokes equations. A good agreement is obtained for small inelasticity and 
system sizes close to the critical value for the shear instability. Finally, Sec. |V] summarizes 
the obtained results and presents a brief discussion. 
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II. HYDRODYNAMIC EQUATIONS 



The balance equations for the number density n{r,t), the flow velocity u{r,t), and the 
granular temperature T(r, t) of a system composed by smooth inelastic hard spheres {d = 3) 
or disks {d = 2) of mass m, diameter a, and coefficient of normal restitution a, have the 
form 

^ + V-(nn) = 0, (1) 
d \ 

— + u-Vju + {mn)-^V-P = 0, (2) 
For a dilute gas and to Navier-Stokes order, the pressure tensor P, the heat flux g, and the 

dnn 

cooling rate ( are given by the constitutive relations [9|, Il7|, [IS] 



P.,(r, t) = p{r, t)5., - + 1^ - l5.,V ■ , (4) 

q(r, t) = -kVT - /iVn, (5) 

C = C^'^ + C^'\ (6) 

Here p = nT is the hydrostatic pressure, rj is the shear viscosity, k the (thermal) heat 
conductivity, and /i a transport coefficient peculiar of granular fluids that is referred to as 
the diffusive heat conductivity. The term (^^^ denotes the contribution to the cooling rate 
of zeroth order in the gradients of the hydrodynamic flelds, while C^^"* denotes the second 
order in the gradients contributions. The latter will be neglected in the following, as it 
is done in most of the calculations, since it is expected to give very small corrections to 
the hydrodynamic equations as compared with the similar terms coming from the pressure 
tensor and the heat flux [o^. The transport coefficients appearing in Eqs. (jl]) and ([5]) can 
be expressed as 

r/(T,a) = r^oW(a), (7) 

K,(T,a) = no(T)K*{a), (8) 

/i(n,T,a) = ^^^/i*(a), (9) 
n 

where tjqIT) and k,q{T) are the low density (Boltzmann) values of the shear viscosity and 
the thermal heat conductivity, respectively. Their explicit expressions as well as those of the 
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dimensionless functions 77*, k*, and /i* are given in Appendix |Xl Finally, the expression of 
the zeroth order cooling rate reads 

C(°Hn,T,a) = ^C(«). (10) 
VoK-i- ) 

The function C*('^) is also given in Appendix Rl 

The hydrodynamic Navier-Stokes equations for the inelastic gas are obtained by substi- 
tuting Eqs. (jll)-(l6]) into Eqs. (IH)-©. They admit a simple solution defined by a constant and 
uniform density uh, a vanishing velocity field uh = 0, and a time dependent temperature 
Tfjit) obeying the equation 

^ + Ci?'(i)r«(i) = o. (11) 

This is the so-called homogeneous cooling state (HCS) ij]. This state is known to be un- 
stable with respect to long wavelength spatial perturbations. More precisely, linear stability 
analysis of the hydrodynamic equations shows that, for large enough systems, the transver- 
sal component of the velocity relative to the square root of the temperature grows in time 
leading to the instability of the system 0,15,01. This shear instability has been confirmed by 
molecular dynamics simulations 0, O] and also by Monte Carlo simulations of the effective 
dynamics associated to the Boltzmann equation 19|]. Of course, after a short time interval, 
the theoretical predictions based on the linearized hydrodynamic equations are not valid any 
longer since neglected nonlinear effects become very important. During the initial stages of 
the development of the instability, density inhomogeneities occur in the system. For dilute 
systems whose linear size is larger (and comparable) to the critical system size for the shear 
instability, the spatial inhomogeneities seem to be dominated by a nonlinear coupling of the 



density with the transversal velocity field 



20| . Although this could indicate that the system 



is going into a clustering regime, in which the particles tend to group together forming very 
high density regions coexisting with very dilute regions, simulation results indicate that the 



20|. 



density actually tends to a quite smooth profile 

Because of continuous cooling due to the inelasticity of collisions, stationary states are not 
possible in freely evolving granular gases (for instance, with periodic boundary conditions). 
Nevertheless, it is still possible that the "final state" reached by the system when the HCS 
is unstable, exhibits some scaling properties so that it can be simply identified at least at 
the hydrodynamic level of description. Here, attention will be focused on systems just above 
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the shear mode instabihty threshold. Well above this threshold, the physical scenario miffht 

npfi 

be quite different [6|, 17|, 18| . 

Consider hydrodynamic ffows verifying the following conditions: i) there are gradients 
in only one direction taken as the X axis, and ii) the hydrodynamic velocity field has the 
form Ui{r, t) = 6i^yu{x, t). Other specifications of the state of the system will be made when 
appropriate. Use of the above conditions into Eqs. ([I])-® shows that the density profile must 
be stationary, n = n{x), and the component Pxx of the pressure tensor must be uniform. 
In the framework of the Navier-Stokes approximation, the latter implies that the pressure 
is also uniform, p = p{t). Moreover, the equations for the velocity and temperature fields 
become 

-(mn)-'-(^-]=0, (12) 

+ TC^°^ = 0. (13) 



dt dx V dx 



9? -2''"" " U 



n/i\ dT 



To solve the above equations, the boundary conditions must be specified. A system of 
particles enclosed in a cubic {d = 3) or square {d = 2) box of side L will be considered. 
Periodic boundary conditions will be assumed at all the boundaries in order to avoid unde- 
sirable wall effects. It will be convenient for the purposes here to introduce a function ip{x) 
by 

n{x) =nH[l + (fi{x)] , (14) 

with uh = N/L'^. The conservation of the number of particles and the periodic boundary 
conditions imply that 

L 

dxip{x) = 0, 

ip{x + L) = ip{x), 0<x<L. (15) 

Moreover, it will be assumed that |v5(a;)| <^ 1, i.e. the spatial density inhomogeneities are 
supposed to be small. It will be shown in the following that this actually reduces the range of 
applicability of the theory to systems near (and above) the threshold of the shear instability. 
Then, the temperature profile in the state being considered is given by 

T{x,t)^9{t)[l-^{x)] (16) 

where ^ 

e{t) = /o dxn{x)T{x,t) ^^^^ 

dxn{x) 
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is the spatial average temperature of the system at time t. It must be reahzed that assuming 
^ 1 and, consistently, retaining terms up to first order in it when solving Eqs. f|T2|) and 
(IT^ . is not the same as linearizing around the HCS, since nothing is in principle assumed 
about the amplitude of the temperature 6{t) or the velocity field u{x,t). This will be 
discussed in detail in the next section. 



III. THE FREE SHEAR STATE 

The hydrodynamic equations fll2l) and (IT^ can be simplified by using dimensionless 
length and time scales defined by 



m 



1/2 



2 [e{t)_ 

and 



X, (18) 



s^\j^ dt'Ht'), (19) 

respectively. Here 

is a characteristic frequency, proportional to the Boltzmann collision frequency ^'^(t), namely 
t^Bit) = (2 + d)v(t)/A. Note that the pre-factor in Eq. ( fTSi) does not depend on time and 
the length scaling is made with a characteristic length which is proportional to the (time- 
independent) average mean free path of the system. In terms of the new scales, Eq. ( lT2l) 
can be approximated by 

ds 2 ^ 

where terms of order u(f have been neglected. The solution of the above equation is a 
superposition of monocromatic waves of the form 

Uk{l, s) = Xfc(0^fc(s) (22) 

with 

Xk{l) = sm{kl + 0fc) (23) 

and 

Ukis) = uJk,oe~^^, (24) 
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where (pk and Ukfi are arbitrary constants. The possible values of k are restricted by the 
periodic boundary conditions which imply 

(25) 

I'M 

q being a positive integer and Im the value of / for x = L, i.e. 



I'M 



m 



1/2 

L. (26) 



2 [e{t)_ 

The form of Eq. ( 12^ indicates that in the limit of large time s, the dynamics of the system 
is governed by the fundamental mode corresponding to the lowest possible value of /c, i.e. 
k = km = '^t^/Im- Therefore, for large enough times, 

= x(/)^(s). (27) 

with %(/) = Xkmi^) and uj{s) = ujk^{s). In the same approximation under consideration, 
|v9(x)| -C 1, the evolution equation for the temperature (fT3|) leads to 

1 de{s) m7]*{a)uj\s) fdxil)y , ^,,^,^0 , .hm , ^ + 2 r,../.^ , ^V(0 n 

eisj^V^J ^"^f' + ^^^^] + ^(^f" («)-M")]^ = o. 

(28) 

Using Eq. (l23l) this becomes 

1 d6'(s) ?7*(Q;)mu;^(s)A;^ 



{l + cos[2(M + 0fe™)]} 



+C(a)[2 + ^(/)] + [«:*(«) - /.*(a)] ^ = 0. (29) 

By requiring the sum of the position dependent terms to cancel it is obtained that 

m = - JtllYit"^ , aco.[2{kJ + 0,J], (30) 
Ad{d + 2)7(fc„, a) 

where 

7(a, fc^) = K (a) - /i (a) - ^^^^-^y^ (31) 

and 

Consistency requires that a be actually independent of s and, moreover, that it be defined 
positive. The former of these conditions together with Eq. fl?I]) leads to 



(32) 



1 de{s) 2 duj{s; ,2 , 



^(s) ds ijj{s) ds 



-klv*{a). (33) 



This equation will be used later on to determine the effective cooling rate of the average 
temperature 9{s) of the system in the free shear state. On the other hand, equating to zero 
the part of Eq. (12^ that is position independent yields, after employing Eq. ([2SD, 

_ 2d[2C{a)-klv*ia)] 

As mentioned above this quantity must be positive. As a consequence, the existence of the 
mathematical solution of the hydrodynamic Navier-Stokes equations under consideration 
and hence of the free shear state is only posible if 

2C 



T]* 



or, equivalently, L > L^, with the critical size given by 

(2 + d)r {d/2) / ^ 



7]^ 

2C 



(35) 



(36) 



This coincides with the condition determining the instability region of the shear mode found 
in the linear stability analysis of the hydrodynamic equations around the HCS Qj. 

Substitution of Eq. flMl) into Eq. fl5U]) provides the explicit form of the steady density 
profile in the free shear state near the threshold of the instability, 

ip{l) = -Acos[2{kml + (l)kJ], 



A 



{d-l)ri*{a) 



^ -1 



(37) 

2(rf + 2)7(fc„,«) ' ' ' ^ ■ ^^^^ 
Therefore, the condition |v2(a^)| ^ 1 formally implies that \L — Lc\/Lc ^ 1. How restrictive 
this condition actually is will be seen in the next section. 

A particularly simple expression for the amplitude of the velocity field is obtained by 
scahng with the average temperature. Equations fl32|) and flMI) yield 



m 



2e{s) 



1/2 



1/2 



(39) 



Therefore, the scaled macroscopic velocity field is time-independent and all its dependence 
on the average density and inelasticity occurs through the critical length L^. A similar 
expression can be obtained for the ratio between the effective cooling rate of the shear state 
Q and the cooling rate of the HCS, C*. The former is identified by writing Eq. fl33l) in the 
form 

de{s) 



ds 



+ 2Cse{s) = 0, 



(40) 



with 



that leads to 



klv*{a) ^^^^ 



(42) 



Of course, in the hmit L — >• Lc, 0{s) — > Th{s) and the law for the HCS given by Eq. (11) is 
recovered. 

The average energy per particle is 

Td 

eT = y«2 + — , (43) 
with the bar denoting spatial average. Using Eqs. ffTTl) . fl22|) . fl23|) . fl32l) . and flM|) it is found: 

For L — > Lc, the equilibrium relation = (iT/2 is recovered as expected. This equation has 
been previously obtained by Wakou et al. 10(], in the context of a Landau-Ginzburg-type 
equation of motion derived from the hydrodynamic equations, under certain restrictions and 
in the limit of nearly elastic collisions (1 — a ^ 1). At this point, it is worth mentioning 
that the assumed periodic boundary conditions do not play an essential role in the theory 
developed here, although they must be compatible with the symmetry of the shear state. 
For instance, it is easily realized that completely equivalent results are obtained if elastic 
walls were used at the boundaries of the system. 

n 

Soto et al., pj| have carried out a nonlinear analysis of the hydrodynamic equations of a 
system of inelastic hard disks close to the instability threshold. When km is slightly smaller 
than kc = 27r/Lc, vorticity modes with wavenumber fc^ exhibit a critical slowing down, so 
that all the other hydrodynamic modes can be considered as enslaved by them. Then, it is 
possible to derive closed equations for the amplitudes of the vorticity modes with k = km 
by using the adiabatic elimination method. This approach can be expected to be formally 
consistent with the one presented here, but quantitative and qualitative discrepancies occur 
when comparing the results from both approaches. Although the authors of fio!] do not 
give almost any detail of their analysis of the nonlinear hydrodynamic equations, we have 
identified a twofold origin of the discrepancies. Firstly, a nonlinear term involving the 
vorticity seems to have been inconsistently neglected. Secondly, the linearization around 



10 



the HCS and the adiabatic method used in [l6j is not equivalent to the approximation 
scheme presented here, in which the velocity field obeys the closed Eq. (|2T1) . 



IV. DIRECT MONTE CARLO SIMULATIONS 

To test the theoretical predictions presented in the previous sections and, in particular, 
the existence itself of the free shear state, we have employed the direct simulation Monte 



Carlo (DSMC) method |21l . |22| | , to numerically generate the dynamics of a system of inelastic 
hard disks. The DSMC method is a many-particle algorithm designed to mimic the effective 
dynamics of the particles of a gas in the low density limit. Therefore, in this limit it is 
expected to lead to the same results as, for instance, molecular dynamics simulations, with 
the advantage of a much larger statistical accuracy. In all the simulations, the initial state 
was homogeneous and isotropic with a Gaussian velocity distribution, and a square box with 
periodic boundary conditions was used. The linear size of the system L was always larger 
than, but close to, the critical one Lc, predicted by Eq. fl5B]) . for the considered value of the 
restitution coefficient a. Then, according to the theory, the system is expected to generate 
the non-linear free shear flow described in Sec. Illlt if it is stable. 

As usual in DSMC simulations, the mass m of the particles will be taken as the unit of 
mass, the average mean free path A = {2\/2nH(y)~^ as the unit of length, and 2T(0) = 26'(0), 
where T(0) is the initial temperature, as the unit of energy. The effect of a collision between 
particles r and s is to instantaneously modify their velocities according to the rule 

1 + a 

Vr ^ V'^ = Vr ■ ""rs)^^ (^5) 

1 + a 

Vs^v'^ = Vs^ —{a-Vrs)a-, (46) 

where Vrs = v^ — Vg is the relative velocity and B is the unit vector pointing from the center 
of particle s to the center of particle r at contact. This corresponds to the scenario in which 
the constitutive relations dl])-© were derived [ol. \l7\. 

One of the technical difficulties when numerically simulating a freely evolving granular 
fluid is the continuous cooling of the system, i.e. the decrease in magnitude of the typical 
velocity of the particles. As a consequence, the numerical inaccuracies become very large 
after some time interval. To avoid this difficulty, a procedure was introduced based on a 
change in the time scale being used. More specifically, a logarithmic time scale is introduced 
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23 



2J]. Then, the dynamics of systems in states where all the time dependence comes 



through the average temperature because of inelastic cooling, can be easily mapped into the 
dynamics around steady states. In this case, the logarithmic time scale is proportional to 
the cumulated number of collisions per particle. 



A. Transient dynamics 

In most of the simulations to be reported, a similar time sequence was observed. The 
system remains homogeneous with no macroscopic velocity field for an initial period of time, 
developing afterwards a state with two vortices, and finally a shear state characterized by 
two counterflows parallel to one of the sides of the system. In this state, there is a coupling 
between the velocity and density fields. The spatial inhomogeneities remain stationary 
and the state looks time independent in the logarithmic time scale used in the numerical 
simulations. The particular location of the transient vortices and their direction, and hence 
of the bands in the shear state, depends on the initial state or, in a statistical sense, on the 
fluctuation taking the system away from the HCS. An example of the described behavior is 
given in Fig. [1], where the scaled velocity field u is plotted at four different times 

for a system with a = 0.95 and L = 39 {L/L^. ^ 1.205). The formation of the vortices, 
their distortion, and the formation of shear bands are clearly identified. The latter remain 
unchanged after their formation, and this was observed in all cases, at least as long as 
the system is close enough to the shear instability. The spontaneous symmetry breaking 
occurring in the final shear state was observed in both perpendicular directions and, to 
compare with the theoretical predictions, the x-axis was always chosen perpendicular to the 
velocity field. More will be commented on this issue in the final section of the paper. 

To get some additional insight about the time evolution of the system before reaching 
the free shear state, in Fig. [2]the ratio between the spatial average temperature 9{t) and the 
HCS temperature at the same time as predicted by the Haff law, Eq. ffTTl) . is shown as a 
function of the average number of collisions per particle r. The system is the same as in Fig. 
[Hand three different simulation trajectories are plotted. Although the time evolution of the 
temperature is not exactly the same in all quite similar trend is observed. There is 

a time interval in which the average temperature is well described by the Haff law (r ^ 250) 
in spite of the system having already well developed vortices (see Fig. [T]). Afterwards, for 
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T=120.5 T=241.2 




10 ' '20 30 ^**nO 20**^'ir*^ 

FIG. 1: Time evolution of the velocity field scaled with the average thermal velocity, tt(m/20)^/^, 
in a system with a = 0.95 and L/Lc ~ 1.205. The indicated times r are measured as the average 
number of collisions per particle. The lengths of the arrows in the plots are proportional to the 
scaled velocities. 

250 ^ r ^ 600, the average temperature decays much slower than in the HCS, so that the 
ratio grows very fast, until it saturates at a given value. The region with the largest growth 
corresponds to the distortion of the velocity vortices. The constant value reached in the 
final shear state indicates that the average temperature 6{t) still obeys a Haff-like law, but 
with a different cooling rate, in qualitative agreement with Eq. (1401) . 

The evolution of the density is illustrated in Fig. [HI also for the system with a = 0.95 
and L/Lc ~ 1.205. Two Fourier components of the relative density field p = n/nn are 
plotted, p2fc„ and Pk„^,km- The former corresponds to the second density mode that is the 
one predicted to survive by the theory, Eq. fl37|) . In the simulations, it appears parallel to 
any of the sides of the system as discussed above and no distinction is made when reporting 
the simulation results. The other component, Pfc„,fc„ is the lowest mode along the diagonal 
of the system. Again, the several curves are different simulation trajectories. In all cases, the 
second transversal mode of the density grows from the noise level to a constant value as the 
shear instability develops. Moreover, the final steady value is the same for all trajectories. 
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200 400 600 800 1000 
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FIG. 2: (Color online) Time evolution of the spatial average temperature normalized by the tem- 
perature of the reference HCS for the same system as in Fig. [TJ Time r is measured as the average 
number of collisions per particle. The several curves correspond to different simulation trajectories. 

On the other hand, although the diagonal mode often also grows in the time interval in 
which the two vortices are losing their shape transforming into the shear state, it decays to 
the noise level once this state is reached. 

B. Shear state results 

Now the results obtained for the properties of the system once in the free shear state will 
be reported, and compared with the theoretical predictions obtained in this paper. Consider 
first the density profile. The simulation data are very well fitted by a cosine function as given 
in Eq. (1371) . This has been checked both by plotting directly the profiles and by computing 
their Fourier components. The amplitude of the measured perturbation An is plotted in Fig. 
Hlas a function of L/L^ for three values of the restitution coefficient, a = 0.97, 0.95, and 0.9. 
The lines are the theoretical predictions given by Eq. fl55]) . the highest one corresponding 
to the smallest value of a and the lowest one to the greatest value of a. It is seen that the 
agreement is quite good when L/Lc is close to 1, the discrepancies increasing as the size of 
the system increases above its critical value. Also, the deviation is larger the smaller the 
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FIG. 3: (Color online) Time evolution of the second transversal mode of the relative density P2fcm 
(solid lines) and the diagonal mode Pkm,km, (dashed lines). The system is the same as in Fig. [TJ 
The several curves correspond to different trajectories, and time r is measured again as the average 
number of collisions per particle. 

coefficient of restitution. In any case, it must be kept in mind that the theory developed 
here is by construction restricted to states with yl„ <^ 1. 

In Fig. [SI the amplitude At of the steady cosine profile for 1 — T{x,t)/9it) is presented 
for the same values of the parameters as considered in Fig. HI According to the theory 
developed here, Eq. (|T6ll . this amplitude should be the same as for the density profile, i.e. 
At = An = a. Nevertheless, the simulation data indicate that this is not the case, except 
for values of L/L^ close to unity. The deviations of the numerical data from the theoretical 
predictions are in opposite directions for At and An, being larger for the former. This 
indicates that the pressure is not actually strictly uniform as predicted by the Navier-Stokes 
equations for dilute granular gases, but exhibits some oscillatory profile. On the other hand, 
the component Pxx of the pressure tensor was found to be uniform, as required by the balance 
equations ([I])- ([3]). Nevertheless, for L/L^ small enough, the agreement between theory and 
simulation can be qualified as satisfactory. Let us mention that if the elastic values of the 
transport coefficients [rj* = n* = 1, fi* = 0) were used, the theoretical prediction for A as 
a function of L/Lc would become independent of a contrary to what is observed in the 
simulations. 



15 



0.3 
0.2 
0.1 


1 1.2 1.4 1.6 

FIG. 4: (Color online) Dimensionless amplitude of the cosine density profile in the free shear state 
as a function of the ratio between the size of the system L and its critical value for the shear 
instability L^- Three values of the restitution coefficient are considered, as indicated. The symbols 
are simulation results and the lines are the predictions of Eq. (j38p . the lowest line corresponding 
to the larger a. 

A quite strong prediction of the theory is provided by Eq. (!39|) . where the amphtude of the 
velocity field scaled with the square root of the average granular temperature is expressed as 
a simple time and a-independent function of the ratio L/Lc- This latter property was seen 
to ve verified in the simulation within the statistical uncertainties. Moreover, the results 
displayed in Fig. O show that the simulation data are in good agreement with the theoretical 
expression, although again systematic deviations are observed as the value of the coefficient 
of restitution decreases and/or the size of the system as compared with its critical value 
increases. 

The simulation results also indicate that the decay of the average temperature in the 
free shear state is governed by a Haff-like law, Eq. (l39l) . In Fig. [71 the ratio between the 
cooling rate of the free shear state, ^5, and the one of the HCS, C*, is plotted, for the same 
systems being considered along this section. The theoretical prediction is provided by Eq. 
(H2|) . implying that the free shear state cools slower than the associated HCS, i.e. with 
the same density and initial temperature. Again a quite satisfactory agreement is found 
between theory and simulations, with the discrepancies exhibiting the same trends as in all 
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FIG. 5: (Color online) The same as in Fig. S] but for the amplitude of the cosine temperature 
profile. 

the previous comparisons. 

V. SUMMARY AND DISCUSSION 

It has been shown that the hydrodjTiamic Navier-Stokes equations for granular gases 
predict the existence of a shear state for freely evolving systems, whose size is slightly larger 
than the critical size for the shear mode instability. Explicit analytical expressions for the 
hydrodynamic fields have been derived. Although it is a non- stationary state due to cooling, 
all the time dependence of the hydrodynamic fields occurs through the average temperature 
and, therefore, it can be eliminated by introducing appropriate dimensionless quantities. 
These theoretical predictions have been found to be in good agreement with the numerical 
results obtained by the DSMC method, for values of the restitution coefficient a close to 
unity. When the system is more inelastic, significant deviations from the theory are observed. 
This was expected, since the free shear state is characterized, as many other non-equilibrium 
states of granular gases, by a strong coupling between gradients and inelasticity. In the 
present case, this is easily identified through the dependence on a of L^. As a consequence, 
a first order gradient expansion like the one leading to the Navier-Stokes equations also 
implies a limitation in the range of values of a for which the theory applies. Moreover, 
it is probably true that the free shear state is inherently non-Newtonian, as it is the case 
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FIG. 6: (Color online) Dimensionless amplitude of the macroscopic flow to = uj{t)/ {26[t)/m) 



1/2 , 



m 



the free shear state as a function of the length L of the system normalized by its critical value Lc- 
The solid line is the theoretical prediction, Eq. (j39|) . while the symbols are simulation results for 
the same systems as in Fig. [H 



for the steady uniform shear state of a granular gas [25^. A possible indication of this is 
the inhomogeneity of the pressure observed in the simulations as implied by the difference 
between the measured values of and At- 

In ref. jcl several possible hydrodynamic scenarios are described for the behavior of a 
freely evolving granular fluid inside its instability region. In that classification, the situation 
considered in this paper is called scenario 4. Here the final state reached by the system has 
been investigated and a quantitative test for the complete scenario has been provided. The 
free shear state seems stable, in the sense that no deviations from it have been observed in 
the simulations. In this context, it is worth mentioning that, in many cases, the two- vortex 
state was clearly identified for a quite large period of time before the system moved to the 
shear state. This may indicate that it is a metastable state with a large escape time. This 
issue clearly deserves more attention. On the other hand, when L becomes much larger than 
Lc (typically, L/Lc ^ 1-7), other hydrodynamic modes become relevant and the free shear 
state is not expected to occur in the system. This has been confirmed by the simulation 
results. 

It is also worth comparing in some detail the approach followed here with the nonlinear 
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FIG. 7: Ratio between the cooling rates of the free shear state and of the HCS as a function of the 
ratio L/Lc- The symbols have the same meaning as in Fig. |3]and the solid line is the theoretical 
prediction, Eq. (jlT|) . 



stability analysis carried out in ref. [IG"]. As already mentioned, the results reported there 
suggest that a term nonlinear in the velocity has been omitted although it is of the same order 
as others that are kept. This term dramatically modifies the results of the stability analysis 
261] . Moreover, in 16| an adiabatic elimination method is used, in which the time derivative 
of the hydrodynamic fields other than the transverse velocity are set equal to zero. In this 
way, these fields can be expressed in terms of u, and when the expressions are inserted into 
the equation for m, a closed equation is obtained for the latter. The approximation followed 
here is different. To lowest order, the transversal velocity field obeys a closed equation by 
itself, Eq. lpTl) . once the function 9{t) has been scaled out by the change of variables. This 
equation differs from the one obtained by the adiabatic elimination method. It is not fully 
clear to us presently the physical origin of this strong discrepancy. 

One relevant question that always arises when using smooth inelastic hard particles to 
model granular gases, is to what extent the results would be modified if more realistic models. 



including for instance rotational degrees of freedom [27j and/or velocity dependent restitu- 
tion coefficients [2], were considered. Although the hydrodynamic Navier-Stokes equations 
including these effects are known in some limiting cases, their analysis is far beyond our 
present reach. Nevertheless, it can be expected that an extension of the energy balance 
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appearing in the simple case discussed here will hold when more dissipation mechanisms 
are included. Then it is our conjecture that a similar shear state but including the new 
rotational effects will also show up. 
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APPENDIX A: NAVIER-STOKES TRANSPORT COEFFICIENTS 



In this Appendix, the explicit expressions of the transport coefficients and the cooling 
rate introduced in Eqs. (Hll-( fTOl) are given. The elastic shear viscosity and heat conductivity 
are 

^0 = (0 vr-V (^T)VV-(-^), (Al) 

1/2 



" 16(rf-l) V2/ 



'{d-l) 



(A2) 



respectively. As usual in the context of granular fluids, the Boltzmann constant has been 
set equal to unity. The factors accounting for the dependence on the restitution coefficient 
are given by 



2d 



Cia) 



K (a) 



[^2*(«)-^C(«)]-Ml + C*(«)], 



/i*(a)=2C(a) 



K*(a) + 



{d-l)c*{a) 
2dC*{a) 



^*(a) = ^(l-a2^ 



2{d-l] 



n -1 



4d 



1 + 
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3c*(a 



z/*(a)-3C(a) 



32 



In the above expressions, 

z/*(a 
1 + a 



(3-3a + 2d)(l + a) 
4d 



1 - 



c*(a) 
64 



d- 1 



d-l 3(d + 8)(l-a) A + 5d - 3(4 - d)a , 

+ — H TTZTT. —C*{a) 
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(A3) 
(A4) 
(A5) 
(A6) 

(A7) 
(A8) 
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32(l-«)(l-2a2) ^^^^ 



9 + 2Ad+ {8d - 41)a + 30a'^{l - a) ' 
It is easily verified that k* and rj* tend to unity when a goes to one, while //* and C* vanish 
in this limit. 
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